model {
for(i in 1:N){
mu_M1[i,1] <- 1  #baseline is the 1st category
mu_M1[i,2] <- exp(alpha1*x[i])
mu_M1[i,3] <- exp(alpha2*x[i])
sum_M[i] <-mu_M1[i,1]+mu_M1[i,2]+mu_M1[i,3]
for (k in 1:3)
{mu_M[i,k] <- mu_M1[i,k]/sum_M[i]}
M[i]~dcat(mu_M[i,])
mu_y[i] <-c*x[i]+beta1*equals(M[i],2)+beta2*equals(M[i],3)
y[i] ~ dnorm(mu_y[i],prec2)

j1[i]~dunif(0.5,100.5)
j2[i]~dunif(0.5,100.5)
j3[i]<- round(j1[i])
j4[i]<- round(j2[i])

mu_M1.2[i,1] <- 1  #baseline is the 1st category
mu_M1.2[i,2] <- exp(alpha1*(x[i]+deltax))
mu_M1.2[i,3] <- exp(alpha2*(x[i]+deltax))
sum_M.2[i] <-mu_M1.2[i,1]+mu_M1.2[i,2]+mu_M1.2[i,3]
for (k in 1:3)
{mu_M.2[i,k] <- mu_M1.2[i,k]/sum_M.2[i]}
mu_y1[i]<- c*(x[i]+deltax)+beta1*mu_M.2[i,2]+beta2*mu_M.2[i,3]
mu_y0[i] <-c*x[i]+beta1*mu_M[i,2]+beta2*mu_M[i,3]
te[i]<-(mu_y1[i]-mu_y0[i])/deltax

mu_M1.3[i,1] <- 1  #baseline is the 1st category
mu_M1.3[i,2] <- exp(alpha1*x[j3[i]])
mu_M1.3[i,3] <- exp(alpha2*x[j3[i]])
sum_M.3[i] <- mu_M1.3[i,1]+mu_M1.3[i,2]+mu_M1.3[i,3]
for (k in 1:3)
{mu_M.3[i,k] <- mu_M1.3[i,k]/sum_M.3[i]}
mu_y2[i]<- c*x[i]+beta1*mu_M.3[i,2]+beta2*mu_M.3[i,3]

mu_M1.4[i,1] <- 1  #baseline is the 1st category
mu_M1.4[i,2] <- exp(alpha1*x[j4[i]])#changed j4 to j3 to reduce variance
mu_M1.4[i,3] <- exp(alpha2*x[j4[i]])#changed j4 to j3 to reduce variance
sum_M.4[i] <- mu_M1.4[i,1]+mu_M1.4[i,2]+mu_M1.4[i,3]
for (k in 1:3)
{mu_M.4[i,k] <- mu_M1.4[i,k]/sum_M.4[i]}
mu_y3[i]<- c*(x[i]+deltax)+beta1*mu_M.4[i,2]+beta2*mu_M.4[i,3]
de[i]<-(mu_y3[i]-mu_y2[i])/deltax
ie[i]<-te[i]-de[i]
}

alpha1 ~ dnorm(0.0,0.01)
alpha2 ~ dnorm(0.0,0.01)
beta1 ~ dnorm(0.0,0.000001)  
beta2 ~ dnorm(0.0,0.000001)  
c ~ dnorm(0.0,0.000001) 
var2 ~ dgamma(1,0.1)
prec2 <-1/var2
}